Evidence of universality for the May-Wigner stability theorem for random networks 
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We consider a random network of nonlinear maps exhibiting a wide range of local dynamics, with the links having normally 
distributed interaction strengths. The stability of such a system is examined in terms of the asymptotic fraction of nodes that 
persist in a non-zero state. Scaling results show that the probability of survival in the steady state agrees remarkably well with 
the May-Wigner stability criterion derived from linear stability arguments. This suggests universality of the complexity-stability 
relation for random networks with respect to arbitrary global dynamics of the system. 



The relation between the structure of a network and 
its dynamical properties has been a problem of long- 
standing importance in many fields, especially in theo- 
retical ecology [1]. A major advance in this area was the 
suggestion by May that the stability of a network can 
be inferred from an analysis of the interactions between 
the network elements [2] . Confining attention only to the 
local stability of an arbitrary equilibrium of the dynam- 
ics, one can ignore explicit dynamics and look at only 
the leading eigenvalues of the linear stability matrix. As- 
suming that the network interactions are random, rigor- 
ous results on the eigenvalue spectra of random matrices 
can be applied [3] . If the stability matrix is comprised of 
elements from a normal distribution with zero mean and 
variance cr^, then the network is almost certainly stable 
if NCa'^ < 1, and unstable otherwise. N is the number 
of nodes in the network and C is the network connectiv- 
ity, i.e., the probability that any two given elements of 
the network are coupled to each other, as reflected in the 
sparsity of the matrix [2,4]. This result is often referred 
to as the May-Wigner stability theorem [5]. 

May's suggestion that increasing network complexity 
leads to decrease in stability was supported by earlier 
numerical simulations [6], but it ran counter to the em- 
pirically established conventional wisdom that biodiver- 
sity promotes ecosystem stability. The original result has 
been criticized on the ground that it is obtained by lin- 
earizing about an assumed equilibrium, and so is inappli- 
cable when either the perturbations from the equilibrium 
are large, or, the dynamics does not settle down to a fixed 
point attractor (e.g., they might undergo periodic oscil- 
lations as in a Lotka-Volterra type system) . The ensuing 
stability vs diversity debate in ecology has resulted in a 
large body of literature attempting to resolve this issue 
one way or another [7] . Although much of the controversy 
may have been due to the methods that different groups 
used to measure complexity and stability [8], and the 
two apparently opposing conclusions have been resolved 
in the specific context of a community assembly model 
[9], the general question of whether network complexity 
is conducive to the long-term persistence of the nodes 
remains unresolved. In addition to ecological networks, 
phenomena where the survival of nodes in a network 



maybe of relevance are power grid breakdown, financial 
market crashes, etc., in short, any system that is sus- 
ceptible to sudden collapse. Further, since the present 
problem is related to the persistence of a trajectory in 
a high-dimensional space with absorbing boundaries, it 
is also of considerable relevance to the general question 
of persistence in non-equilibrium systems that has seen 
a huge spurt of interest recently [10]. 

In this paper, we report results on the role that net- 
work complexity plays on global stability (in contrast to 
local stability) of a network, by looking at the persis- 
tence of individual nodes in a network of randomly cou- 
pled nonlinear maps undergoing a wide range of local dy- 
namics. We observe that the results of the May-Wigner 
theorem seem to be valid universally, namely, increasing 
the number of interactions per node or increasing inter- 
action strength will give rise to increased likelihood of 
extinction. This evidence of universality (in the sense of 
being independent of the local dynamics at the nodes) 
has bearing on network problems in general [11,12], as 
it addresses an issue which arises in many different con- 
texts, namely: what is the significance of local dynamics 
on network stability, especially in situations where the 
dynamics can he widely varying. 

Previous work on including explicit dynamics in net- 
work models mostly involved generalized Lotka-Volterra 
type ordinary difi^erential equations (ODEs) [13]. How- 
ever, in the absence of interaction between the nodes, the 
local dynamics in such a system is trivial. In contrast, 
considering randomly coupled maps as a model for the 
dynamical network allows us to consider very general lo- 
cal dynamics, including chaos. In the specific context of 
ecological networks, this is a reasonable assumption for 
the population dynamics of individual species. In addi- 
tion, the use of coupled maps allow us to work with much 
larger networks, compared to models incorporating re- 
alistic consumer-resource configurations used to analyze 
simple communities with very few species, whose results 
are difficult to scale to larger ecosystems [14]. 

Our model has N dynamical elements in a network 
with random nonlocal connectivity, for instance repre- 
senting an ecological network of N interacting species. 
Each node i{= 1 . . . N) is associated with a continuous 
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state variable Xn{i) which represents the relative pop- 
ulation density of the ith species at time n. The in- 
teraction between two species is represented by Lotka- 
Volterra type relation, with the sign of the coupling co- 
efficient Jij determining either a predator-prey relation 
{Jij > 0,Jji < 0), competition {Jij,Jji < 0) or mutual- 
ism {Jij,Jji > 0). The time-evolution of the system is 
given by 

Xn+l{i) = f[Xn{i){'^ + ^jJijXnU)}], (1) 

where / represents the local on-site dynamics. For the 
results shown in this paper we have chosen / to be the 
exponential map, 

/(x) = a;e''(^-^\ if a; > 0; = 0, otherwise. (2) 

r being the nonlinearity parameter leading from periodic 
behavior to chaos [15]. This is a much more realistic 
model of population dynamics than the logistic map, and 
in contrast to the latter, is defined over the semi-infinite 
interval [0, oo] rather than a finite, bounded interval. Our 
results also hold for other models of population dynamics 
such as the Bellows map, f{x) = rx/ (1 -|- x^) [16]. These 
maps have the property that they do not go extinct in the 
absence of coupling, as we are interested not in intrinsic 
instability of the species, but rather in the instability 
induced by network interactions. 

The connectivity matrix J = {Jij} is a sparse matrix 
with probability 1 — C that an element is zero. The di- 
agonal entries Ju = indicate that in the absence of 
interaction with other species, the exponential map (2) 
completely determines the population dynamics of each 
species. The non-zero entries in the matrix are chosen 
from a normal distribution with mean and variance . 
Note that we have also used uniform distribution over 
the interval [— ct, a] without any qualitative changes in 
the results. The results reported below are for parallel 
updating; similar results hold for random sequential up- 
dating. Also, our results hold for interaction couplings 
other than the one used above. For example, the follow- 
ing type of coupling: 

Xn+l{i) = f[Xn{i)] + '^jJijf[Xn{i)]f[Xn{j)], 

gives results similar to that reported in this paper. 

The linear stability criteria for random networks pro- 
vides a relation between the parameters N, C and a. 
However, since we are considering explicit local dynam- 
ics, we have an additional parameter, r. In our work, 
instead of looking at linear stability, we shall consider 
persistence, i.e., the probability that a site has a non- 
zero value of the measure of stability of the system. 
Although some early work on survival and extinction of 
species in a coupled network were done in restricted con- 
texts of exclusively competitive [17] or cooperative in- 
teractions [18], no systematic study has been previously 
made on whether the May criterion is valid in the pres- 
ence of local dynamics, incorporating all kinds of inter- 
actions between species. 



Initially, all the N species have population values ran- 
domly distributed about x = 1. Immediately after start- 
ing the simulation the number of persistent species (i.e., 
with a; > 0) decreases rapidly, but eventually attains a 
steady state value which is a function of the system pa- 
rameters. Note that, if a; < for any species, it is re- 
moved from the system and subsequently plays no fur- 
ther role. After a series of such extinctions, the effec- 
tive number of interacting species decreases and, conse- 
quently, the intensity of such extinction-inducing fluctua- 
tions is also reduced. We have continued the simulations 
for up to 10^ iterations, when the probability of further 
extinctions was found to become extremely small. We 
then look at the fraction of species which survive as a 
function of the model parameters (Fig. 1). The results 
qualitatively agree with the May criterion for stability, 
in that, increasing complexity (in terms of size, connec- 
tivity and interaction strength of the network) decreases 
stability, with a larger proportion of species liable to get 
extinct. Note that the May criterion was derived on the 
basis of local stability, whereas here we are considering 
the species persistence, a measure of global stability. 

Fig. 1(a) shows the ratio of persistent species Npers 
with respect to the initial number of species N. This 
ratio Npers/N appears to vary as 1/N for large N. This 
indicates that the number of surviving species is indepen- 
dent of N. Agreement with Wigner-May stability results 
is also seen for the 1/C variation of surviving fraction 
with connectivity (Fig. 1(b)). Fig. 1(c) shows that the 
fraction of survivors depend on the interaction strength 
parameter cr as 1 /a^ where the exponent z is an increas- 
ing function of the connectivity C. This dependence is 
expected because, if C is decreased while keeping N fixed, 
the effective number of other species that a species inter- 
acts with, is decreased. In the limit C ^ 0, every species 
is independent of all other species, and will persist with 
probability 1. Finally, we display the survival fraction 
against the nonlinearity parameter r of the local map. It 
is clearly evident that one obtains a smooth monotonic 
variation of the survival fraction with respect to r (Fig. 
1(d)). This a priori may seem surprising, since the local 
map has a significant range of diverse dynamics including 
windows of periodic and chaotic behavior and this is not 
reflected at all in the figure. 

To understand these results, we analyze the probabil- 
ity of survival of any species in the steady state. A species 
i will become extinct if its population Xi becomes nega- 
tive at a particular time. By looking at the equations de- 
scribing the system, one notes that this is only possible if 
T,jJijXj < — 1. Therefore, the probability of survival of a 
species is essentially equivalent to P(I]j JijXj > —1). The 
distribution of P{Y,jJijXj) has a power law distribution 
about its peak at zero, and Gaussian tails. We now scale 
this distribution with respect to the different network pa- 
rameters, as scaling in non-equilibrium phenomena is the 
most sensitive and stringent test of universality. 

Fig. 2 shows the scaling of P{Y,jJijXj) with the con- 
nectivity C which goes as ~ C~^gc{C~^'EjJijXj) where 
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Qc is the scaling function independent of C, implying 
CP{T,jJijXj > —1) ~ constant. Therefore, the proba- 
bility of survival varies as ~ 1/C, in exact agreement 
with the results obtained from linear stability analysis. 
The exponent /? = 0.2 ±0.02 for a wide range of values of 
a and r. Similar agreement is seen for the variation of the 
probability of survival with a (Fig. 3). The scaling data 
show that P{TijJijXj) ^ a~'^ga-ia'~°'^jJijXj), where 
is the scaling function, so that the survival probability 
varies as The exponent a varies in the range 

0.1-0.2, decreasing with r and with C. 

The variation with the map nonlincarity parameter 
however has no analog in the previous work on random 
networks. We observe that the relevant parameter is the 
image of the critical point of the map, rather than r it- 
self. This point a;™°^ = e^^'^^'/r gives a measure of the 
width of the chaotic attractor [19]. Since this increases 
the interval over which the probability of (TijJijXj) is 
observed, we have normalized the argument of the scal- 
ing function by dividing it by x^"''^. Fig. 4 shows the 
scaling of P(Ej Jyx,) - (x™«^)-T54(x™»-)-i5],- Ji,-Xj], 
where is the scaling function. Therefore, the proba- 
bility of survival varies as (a;™"^ )"''', with the exponent 
7 = 3.1 ± 0.1 for a wide range of values of C and a. Inter- 
estingly, when the local dynamics is given by the Bellows 
map, we again obtain 7 ~ 3. 

The above scaling results show that the complexity- 
stability relations obtained by May hold true not only 
qualitatively, but also quantitatively, when we introduce 
explicit local dynamics of the network elements. The 
exact nonlinearity of the map, as would be reflected in, 
e.g., the Lyapunov exponent, does not enter any of the 
results, which suggests that these relations arc^ univer- 
sal and independent of details of the local dynamics. In 
addition, the results remain valid even if the local nonlin- 
earity parameter r for all the N maps is not a constant, 
but varies according to a uniform random distribution 
between r — 2 and r = 4. 

The power spectra of quantities such as the total sys- 
tem population, X^^j^ Xi (which can be identified with 
"biomass" in the ecological context), has a low frequency 
scaling given by : S{f) ~ /~" with 1 < a < 2. In 
addition, the distribution of populations P{x) is a clear 
power law: P{x) ^ x~^, with <j) 1^ 1 for sufficiently high 
r [Fig. 4 (inset)] [20]. 

In summary, our work addresses one of the strong crit- 
icisms against the wider applicability of the May-Wigner 
results, namely their assumption of an equilibrium. Here 
we have a range of dynamics at the local level and cer- 
tainly no dynamical equilibrium at the global level, as 
populations are always fluctuating. Rather we have a 
non-equilibrium steady state where the survival fraction 
attains stationarity. The stability of our dynamically 
more complex network however still obeys the May cri- 
terion, and increasing complexity (in terms of size, con- 
nectivity and interaction strength of the network) leads 
to greater instability, resulting in a larger proportion of 



species becoming extinct [21] . Scaling results of the prob- 
ability distribution of the interaction term in the sta- 
tionary state indicate that the stability of the network 
varies as ~ MCa'^ ' ^'^^^ much in agreement with the May- 
Wigner results. We also find that the stability of the net- 
work scales with the nonlinearity parameter of the local 
maps in a smooth monotonic fashion, with the relevant 
scaling variable being the maximum value that x can 
take (which depends monotonically on the nonlinearity). 
These observations hold for networks with widely vary- 
ing local dynamics as well as for different updating and 
coupling schemes, underscoring a remarkable universality 
and increasing the scope of relevance of the May-Wigner 
stability theorem. 
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FIG. 1. The fraction of persistent nodes plotted against 
the model parameters: (a) the initial number of nodes, A'^ 
(cr = 0.1; o: C = 0.1, r = 2, □: C = 0.1, r = 4, A: 
C = 1, r = 2, v: C = 1, r = 4; (b) connectivity, C {N = 100, 
a = 0.15; o: r = 2, □: r = 3, A: r = 4); (c) standard 
deviation, a {N ^ 100, r = 4; A: C = 0.1, □: C = 0.25, 
\j: C = 0.5, o: C = 1); (d) the nonlinearity parameter, r 
{N = 100, a = 0.1; o: C = 0.1, □: C = 0.5, A: C = 0.9). 
The data is obtained after 10* iterations and averaged over 
5000 realizations. 
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FIG. 3. The scaling of T,jJijXj with cr, the standard devia- 
tion of normal distribution from which the connection weights 
are chosen [N = 100, C = 1 and r = 4). The data is obtained 
after 10* iterations and averaged over 5000 realizations. 
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FIG. 4. The scaling of T.jJijXj with the width of the at- 
tractor a;™"'' for N = 100, C = 0.5 and cr = 0.1. The inset 
shows the power-law scaling behavior of the probability dis- 
tribution of populations x {N = 100, C — 1 and cr = 0.1). 
The data is obtained after 10* iterations and averaged over 
5000 realizations. 



FIG. 2. The scaling of 'EjJijXj with connectivity C for 
= 100, cr = 0.1 and r = 4. The data is obtained after 
10* iterations and averaged over 5000 realizations. 
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